//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics 
//
//  License:		 BSD License 
//					 Kratos default license: kratos/license.txt
//
//  Main authors:    Author Julio Marti
//

// Project includes 
#include "includes/define.h"
#include "custom_conditions/rad_face3D.h"
#include "utilities/math_utils.h"
#include "radiation_application.h"
#include "includes/kratos_flags.h"
#include "includes/deprecated_variables.h"


namespace Kratos
{
  //************************************************************************************
  //************************************************************************************
  RadFace3D::RadFace3D(IndexType NewId, GeometryType::Pointer pGeometry)
    : Condition(NewId, pGeometry)
  {		
    //DO NOT ADD DOFS HERE!!!
  }
  
  //************************************************************************************
  //************************************************************************************
  RadFace3D::RadFace3D(IndexType NewId, GeometryType::Pointer pGeometry,  PropertiesType::Pointer pProperties)
    : Condition(NewId, pGeometry, pProperties)
  {
  }
  
  Condition::Pointer RadFace3D::Create(IndexType NewId, NodesArrayType const& ThisNodes,  PropertiesType::Pointer pProperties) const
  {
    return Condition::Pointer(new RadFace3D(NewId, GetGeometry().Create(ThisNodes), pProperties));
  }
  
  RadFace3D::~RadFace3D()
  {
  }
  
  //************************************************************************************
  //************************************************************************************
  void RadFace3D::CalculateRightHandSide(VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
  {
    //calculation flags
    bool CalculateStiffnessMatrixFlag = false;
    bool CalculateResidualVectorFlag = true;
    MatrixType temp = Matrix();
    
    CalculateAll(temp, rRightHandSideVector, rCurrentProcessInfo, CalculateStiffnessMatrixFlag, CalculateResidualVectorFlag);
  }
  
  //************************************************************************************
  //************************************************************************************
  void RadFace3D::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
  {
    //calculation flags
    bool CalculateStiffnessMatrixFlag = true;
    bool CalculateResidualVectorFlag = true;
    
    CalculateAll(rLeftHandSideMatrix, rRightHandSideVector, rCurrentProcessInfo, CalculateStiffnessMatrixFlag, CalculateResidualVectorFlag);
  }
  
  //************************************************************************************
  //************************************************************************************
  void RadFace3D::CalculateAll(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, 
			       ProcessInfo& rCurrentProcessInfo,
			       bool CalculateStiffnessMatrixFlag,
			       bool CalculateResidualVectorFlag)
  {
    KRATOS_TRY
      
      unsigned int number_of_nodes = GetGeometry().size();
    
    //resizing as needed the LHS
    unsigned int MatSize=number_of_nodes;
    
    //calculate area
    double area = GetGeometry().Area();
    
    const double& ambient_temperature = GetProperties()[AMBIENT_TEMPERATURE];
    double StefenBoltzmann = 5.67e-8;
    double emissivity = GetProperties()[EMISSIVITY];
    double convection_coefficient = GetProperties()[CONVECTION_COEFFICIENT];
    
    const double& T0 = GetGeometry()[0].FastGetSolutionStepValue(TEMPERATURE);
    const double& T1 = GetGeometry()[1].FastGetSolutionStepValue(TEMPERATURE);
    const double& T2 = GetGeometry()[2].FastGetSolutionStepValue(TEMPERATURE);
    
    const double& q0 = GetGeometry()[0].FastGetSolutionStepValue(FACE_HEAT_FLUX);
    const double& q1 = GetGeometry()[1].FastGetSolutionStepValue(FACE_HEAT_FLUX);
    const double& q2 = GetGeometry()[2].FastGetSolutionStepValue(FACE_HEAT_FLUX);
    
    if (CalculateStiffnessMatrixFlag == true) //calculation of the matrix is required
      {
	if(rLeftHandSideMatrix.size1() != MatSize )
	  rLeftHandSideMatrix.resize(MatSize,MatSize,false);
	noalias(rLeftHandSideMatrix) = ZeroMatrix(MatSize,MatSize);
	rLeftHandSideMatrix(0,0) = ( convection_coefficient + emissivity*StefenBoltzmann*4.0*pow(T0,3)  )* 0.333333333333 * area;
	rLeftHandSideMatrix(1,1) = ( convection_coefficient + emissivity*StefenBoltzmann*4.0*pow(T1,3)  )* 0.333333333333 * area; 
	rLeftHandSideMatrix(2,2) = ( convection_coefficient + emissivity*StefenBoltzmann*4.0*pow(T2,3)  )* 0.333333333333 * area; 
      }

    //resizing as needed the RHS
    double aux = pow(ambient_temperature,4);
    if (CalculateResidualVectorFlag == true) //calculation of the matrix is required
      {
	if(rRightHandSideVector.size() != MatSize )
	  rRightHandSideVector.resize(MatSize,false);
	
	rRightHandSideVector[0] =  emissivity*q0 - emissivity*StefenBoltzmann*(pow(T0,4) - aux)   
	  -  convection_coefficient * ( T0 - ambient_temperature);
	
	rRightHandSideVector[1] =  emissivity*q1  - emissivity*StefenBoltzmann*(pow(T1,4) - aux)  
	  -  convection_coefficient * ( T1 - ambient_temperature);
	
	rRightHandSideVector[2] =  emissivity*q2  - emissivity*StefenBoltzmann*(pow(T2,4) - aux)  
	  -  convection_coefficient * ( T2 - ambient_temperature);
	
	rRightHandSideVector *= 0.3333333333333*area;
	
      }
    //KRATOS_WATCH(area);
    //KRATOS_WATCH(rLeftHandSideMatrix);
    //KRATOS_WATCH(rRightHandSideVector);
    KRATOS_CATCH("")
      }
  
  
  
  //************************************************************************************
  //************************************************************************************
  void RadFace3D::EquationIdVector(EquationIdVectorType& rResult, ProcessInfo& CurrentProcessInfo)
  {
    unsigned int number_of_nodes = GetGeometry().PointsNumber();
    if(rResult.size() != number_of_nodes)
      rResult.resize(number_of_nodes,false);
    for (unsigned int i=0;i<number_of_nodes;i++)
      {
	rResult[i] = (GetGeometry()[i].GetDof(INCIDENT_RADIATION_FUNCTION)).EquationId();
      }
  }
  
  //************************************************************************************
  //************************************************************************************
  void RadFace3D::GetDofList(DofsVectorType& ConditionalDofList,ProcessInfo& CurrentProcessInfo)
  {
    ConditionalDofList.resize(GetGeometry().size());
    for (unsigned int i=0;i<GetGeometry().size();i++)
      {
	ConditionalDofList[i] = (GetGeometry()[i].pGetDof(INCIDENT_RADIATION_FUNCTION));
      }
  }
  
} // Namespace Kratos


